function res = hydro(GRID_SIZE, L, R_SIZE, Dt2, tfinal, GAUSS_AMP, GAUSS_SIZE)
  % switches
  HORIZON_ON = false;
  METRIC_ON = true;


  % dirty tricks to make arguments global
  global Dt DX;
  Dt = Dt2;

  % First order hydro
  %  eliminating the u-dot-dot terms 
  %  fsolve search method

  % prepare figures for plotting
  close all;
  fighydro  = figure; 
  figmetric = figure;
  fighorizon = figure;

  % load color maps
  load('cm', 'mycmap');
  CMAP_OFFSET = 512;
  set(fighydro, 'Colormap', mycmap);
  set(figmetric, 'Colormap', mycmap);

  % compute constants & the r dimension
  DX = L / GRID_SIZE;
  TIME_STEPS = round(tfinal / Dt);
  VECTOR_SUBSAMPLING = 16;
  [gridx, gridy] = meshgrid([1:VECTOR_SUBSAMPLING:GRID_SIZE], [1:VECTOR_SUBSAMPLING:GRID_SIZE]);

  % reserve space for the Gdd computation
  ressize = TIME_STEPS*GRID_SIZE*GRID_SIZE*R_SIZE*4*4*4;
  %disp(sprintf('Data size   : %.2fMB', ressize/1024/1024));
  %res = zeros(TIME_STEPS, GRID_SIZE, GRID_SIZE, R_SIZE, 4, 4);
  res = zeros(TIME_STEPS);  % integrated Ricci scalar

  % R --- coordinate values
  global R;
  rmag = 5;
  R = rmag * tan([0:R_SIZE-1]*pi/2/R_SIZE+0.01);
  R

  rticks = transpose(linspace(0,R_SIZE-1,6));
  rticklabels = num2str( rmag * tan(rticks*pi/2/R_SIZE), '%.3g');


  % summary
  disp(sprintf('Grid Size   : %d', GRID_SIZE));
  disp(sprintf('Grid Length : %f', L));
  disp(sprintf('Cell Length : %f', DX));
  disp(sprintf('Target Time : %f', tfinal));
  disp(sprintf('Time step   : %f', Dt));
  disp(sprintf('Total steps : %d', TIME_STEPS));
  pause(1);


  % video stuff
  VIDEO = false;
  EVERY_FRAME=true;
  FPS = 30;
  SECONDS_PER_LIGHTYEAR = 5;
  if (VIDEO)
    vt = VideoWriter('metric.avi', 'Motion JPEG AVI');
    vt.Quality = 100;
    vt.FrameRate = 30;
    open(vt);
  end




  global A;
  A = sqrt(3)*pi/18 - log(3)/2;

  T  = zeros(GRID_SIZE, GRID_SIZE);
  ux = zeros(GRID_SIZE, GRID_SIZE);
  uy = zeros(GRID_SIZE, GRID_SIZE);
  ut =  ones(GRID_SIZE, GRID_SIZE);
  P1 = zeros(GRID_SIZE, GRID_SIZE);
  P2 = zeros(GRID_SIZE, GRID_SIZE);

  
  Tdt  = zeros(GRID_SIZE, GRID_SIZE); 
  uxdt = zeros(GRID_SIZE, GRID_SIZE); 
  uydt = zeros(GRID_SIZE, GRID_SIZE); 
  utdt = zeros(GRID_SIZE, GRID_SIZE);
  P1dt = zeros(GRID_SIZE, GRID_SIZE);
  P2dt = zeros(GRID_SIZE, GRID_SIZE);

  Tdx = zeros(GRID_SIZE, GRID_SIZE); 
  Tdy = zeros(GRID_SIZE, GRID_SIZE); 
  uxdx = zeros(GRID_SIZE, GRID_SIZE); 
  uxdy = zeros(GRID_SIZE, GRID_SIZE); 
  uydx = zeros(GRID_SIZE, GRID_SIZE); 
  uydy = zeros(GRID_SIZE, GRID_SIZE); 
  P1dx = zeros(GRID_SIZE, GRID_SIZE); 
  P1dy = zeros(GRID_SIZE, GRID_SIZE); 
  P2dx = zeros(GRID_SIZE, GRID_SIZE); 
  P2dy = zeros(GRID_SIZE, GRID_SIZE); 

  u = zeros(GRID_SIZE, GRID_SIZE, 3, 4);

  utdtdt = zeros(GRID_SIZE, GRID_SIZE);
  utdtdx = zeros(GRID_SIZE, GRID_SIZE);
  utdtdy = zeros(GRID_SIZE, GRID_SIZE);
  utdxdx = zeros(GRID_SIZE, GRID_SIZE);
  utdxdy = zeros(GRID_SIZE, GRID_SIZE);
  utdydy = zeros(GRID_SIZE, GRID_SIZE);

  uxdtdt = zeros(GRID_SIZE, GRID_SIZE);
  uxdtdx = zeros(GRID_SIZE, GRID_SIZE);
  uxdtdy = zeros(GRID_SIZE, GRID_SIZE);
  uxdxdx = zeros(GRID_SIZE, GRID_SIZE);
  uxdxdy = zeros(GRID_SIZE, GRID_SIZE);
  uxdydy = zeros(GRID_SIZE, GRID_SIZE);

  uydtdt = zeros(GRID_SIZE, GRID_SIZE);
  uydtdx = zeros(GRID_SIZE, GRID_SIZE);
  uydtdy = zeros(GRID_SIZE, GRID_SIZE);
  uydxdx = zeros(GRID_SIZE, GRID_SIZE);
  uydxdy = zeros(GRID_SIZE, GRID_SIZE);
  uydydy = zeros(GRID_SIZE, GRID_SIZE);

  global g gd g_old gi;
  global detg gr;
  global gam gamd gam_old;
  global GUd Gdd;
  g       = zeros(GRID_SIZE, GRID_SIZE, R_SIZE, 4, 4);
  g_old   = zeros(GRID_SIZE, GRID_SIZE, R_SIZE, 4, 4);
  detg    = zeros(GRID_SIZE, GRID_SIZE, R_SIZE);
  gr      = zeros(GRID_SIZE, GRID_SIZE, R_SIZE);
  gam     = zeros(GRID_SIZE, GRID_SIZE, R_SIZE, 4, 4, 4);
  gam_old = zeros(GRID_SIZE, GRID_SIZE, R_SIZE, 4, 4, 4);
  GUd     = zeros(GRID_SIZE, GRID_SIZE, R_SIZE, 4, 4);
  Gdd     = zeros(GRID_SIZE, GRID_SIZE, R_SIZE, 4, 4);





  sumT = 0;
  % initial conditions for vx vy come here
  % vx
  % vy
  % T
  TEMP_OFFSET =50;
  cx = 0.5 * GRID_SIZE + 0.5;
  cy = 0.5 * GRID_SIZE + 0.5;
  gaussian_size = GAUSS_SIZE * L;
  sig2 = (gaussian_size/DX).^2 ;
  for x = 1:GRID_SIZE
    for y = 1:GRID_SIZE
      T(x, y) =  GAUSS_AMP .* exp(-((x-cx)^2 + (y-cy)^2)/sig2);
      %T(x,y) = 0.2*sin(2*pi*x/GRID_SIZE);
    end
  end
  T = T + TEMP_OFFSET;



  % ****************************************************************
  %  compute derivative kernels
  % ****************************************************************
  % can upgrade on this -- trefethen pg 6

  global d1 d2;
  d1 = zeros(GRID_SIZE, GRID_SIZE);
  d2 = zeros(GRID_SIZE, GRID_SIZE);
  for x = 1:GRID_SIZE
    for y=1:GRID_SIZE
      s1 = 0;
      s2 = 0;
      for k = -GRID_SIZE/2 : GRID_SIZE/2
        s1 = s1 + 2. * pi * i / GRID_SIZE / GRID_SIZE * k * exp(2 * pi * i * k * (x - y) / GRID_SIZE);
        s2 = s2 + (2. * pi * i / GRID_SIZE * k)^2 / GRID_SIZE * exp(2 * pi * i * k * (x - y) / GRID_SIZE);
      end
      d1(x,y) = real(s1)/DX;
      d2(x,y) = real(s2)/DX/DX;
    end
  end
  h = 2*pi/GRID_SIZE;
  col1 = [0 .5*(-1).^(1:GRID_SIZE-1).*cot((1:GRID_SIZE-1)*h/2)]  * h/DX;
  max(max(d1 - toeplitz(col1, col1([1 GRID_SIZE:-1:2]))))
  col2 = [(-pi^2/(3*h*h)-1/6) .5*(-1).^(0:GRID_SIZE-2).*(csc((1:GRID_SIZE-1)*h/2).^2)]  * h*h/(DX*DX);
  % d2 = toeplitz(col2);
  


  Tdx = dx(T);
  Tdy = dy(T);

  % initial conditions for u = 0
  %P1 = (-1/4).*T.^(-1).*((-1).*Tdx.^2+(-1).*A.*Tdx.^2+Tdy.^2+A.*Tdy.^2);
  %P2 = (-1/2).*T.^(-1).*((-1).*Tdx.*Tdy+(-1).*A.*Tdx.*Tdy);

  %P2 = (1/4).*(1+A).*(dy(Tdy)+dx(Tdx));
  %P1 = (1/4).*(1+A).*(dy(Tdy)+dx(Tdx));

  P1dx = dx(P1);
  P1dy = dy(P1);
  P2dx = dx(P2);
  P2dy = dy(P2);

    t = 1;
    phydro();
    pause(0.003);
  


  for t = 1:TIME_STEPS
    tic;

    % ids of the past values of u
    uid = rem(t+3:-1:t,4) + 1;

    % solving non-linear equation of type:
    % L(u, T, u-dot, T-dot) + R(u, T) == 0  -- 3 components
    % each of L, R has 3 components 

    [Tdt uxdt uydt P1dt P2dt] = h2eom(T, ux, uy, P1, P2, Tdx, Tdy, uxdx, uxdy, uydx, uydy, P1dx, P1dy, P2dx, P2dy);


    T  = T  + Dt * Tdt;
    ux = ux + Dt * uxdt;
    uy = uy + Dt * uydt;
    P1 = P1 + Dt * P1dt;
    P2 = P2 + Dt * P2dt;

    ut = sqrt(1+ux.^2+uy.^2);
    u(:,:,1,uid(1)) = ut;
    u(:,:,2,uid(1)) = ux;
    u(:,:,3,uid(1)) = uy;

    % force the time derivatives to zero for the very first time step
    if t == 1
      for k=1:4
        u(:,:,1,k) = ut;
        u(:,:,2,k) = ux;
        u(:,:,3,k) = uy;
      end
    end

    utdt = (u(:,:,1,uid(1)) - u(:,:,1,uid(2)))./Dt;



    Tdx = dx(T);
    Tdy = dy(T);
    utdx = dx(ut);
    utdy = dy(ut);
    uxdx = dx(ux);
    uxdy = dy(ux);
    uydx = dx(uy);
    uydy = dy(uy);
    P1dx = dx(P1);
    P1dy = dy(P1);
    P2dx = dx(P2);
    P2dy = dy(P2);

    


    if (METRIC_ON)    
      % compute all second order derivatives

      % back differencing for dtdt
      utdtdt = (2.*u(:,:,1,uid(1)) - 5.*u(:,:,1,uid(2)) + 4.*u(:,:,1,uid(3)) - 1.*u(:,:,1,uid(4)))./Dt./Dt;
      uxdtdt = (2.*u(:,:,2,uid(1)) - 5.*u(:,:,2,uid(2)) + 4.*u(:,:,2,uid(3)) - 1.*u(:,:,2,uid(4)))./Dt./Dt;
      uydtdt = (2.*u(:,:,3,uid(1)) - 5.*u(:,:,3,uid(2)) + 4.*u(:,:,3,uid(3)) - 1.*u(:,:,3,uid(4)))./Dt./Dt;

      % spectral for all others
      utdtdx = dx(utdt);
      uxdtdx = dx(uxdt);
      uydtdx = dx(uydt);

      utdtdy = dy(utdt);
      uxdtdy = dy(uxdt);
      uydtdy = dy(uydt);

      utdxdx = dx(utdx);
      uxdxdx = dx(uxdx);
      uydxdx = dx(uydx);

      utdxdy = dy(utdx);
      uxdxdy = dy(uxdx);
      uydxdy = dy(uydx);

      utdydy = dy(utdy);
      uxdydy = dy(uxdy);
      uydydy = dy(uydy);
       
      % call the metric code
      metric(T, ut, ux, uy, utdt, uxdt, uydt, utdx, utdy, uxdx, uxdy, uydx, uydy, ...
          utdtdt, utdtdx, utdtdy, utdxdx, utdxdy, utdydy, ...
          uxdtdt, uxdtdx, uxdtdy, uxdxdx, uxdxdy, uxdydy, ...
          uydtdt, uydtdx, uydtdy, uydxdx, uydxdy, uydydy);

      % store the einstein tensor
      % res(t, :, :, :, :, :) = Gdd(:,:,:, :, :);
      
      % ************************************************
      %   check if satisfy Einstein Field Equations
      % ************************************************

      G = zeros(GRID_SIZE ,GRID_SIZE, R_SIZE);
      % in 3+1 dim, Gdd = Rdd - (1/2R + 3) gdd == 0, so that R - (1/2R + 3) 4 == 0, so that R = -6
      % Thus G = -R - 12

      G = gr - 12;

      % display the result
      gr = G;
      
      % integrate |G| with L2 norm;
      %  square G
      G = G.*G;
      % sum along x and y directions
      Gixy = zeros(1,R_SIZE);
      for ii=1:4
        for jj=1:4
          Gixy = Gixy + reshape(G(ii,jj,:), [1,R_SIZE]);
        end
      end
      Gixy = Gixy .* DX .* DX;;

      % integrate along the r direction
      Gixyr = 0;
      for ir=3:R_SIZE-1
        Gixyr = Gixyr + 0.5 .* (Gixy(ir) + Gixy(ir+1)) .* ( R(ir+1) - R(ir) );
      end

      res(t) = Gixyr;
      Gixyr




    end

    if (HORIZON_ON)
      horizon = zeros(GRID_SIZE, GRID_SIZE);
      for xx=1:GRID_SIZE
        for yy=1:GRID_SIZE
          horizon(xx,yy) = fzero(@fgtt, 5, optimset('TolX', 1e-2));
        end
      end
    end
    
    %if rem(t, 10) == 0
      phydro();
      if (METRIC_ON)
        pmetric2();
      end
      if (HORIZON_ON)
        phorizon();
      end
      refresh();
    %end
    disp(sprintf('time step: %04d, total time: %f,  step realtime: %8.4fs', t, t*Dt, toc));


  end














  function res = fgtt(rval)
    f = (1-(T(xx,yy)./rval).^3);
    res = -rval.^2.*f.*ut(xx,yy).*ut(xx,yy) + rval.^2.*(-1 + ut(xx,yy).*ut(xx,yy));
    delu = utdt(xx,yy)+uxdx(xx,yy)+uydy(xx,yy);
    res = res + rval.*delu.*ut(xx,yy).*ut(xx,yy);
    res = res -2.*rval.*ut(xx,yy).*(ut(xx,yy).*utdt(xx,yy)+ux(xx,yy).*utdx(xx,yy)+uy(xx,yy).*utdy(xx,yy));

    br = rval./T(xx,yy);
    F = -sqrt(3)./3.*atan(sqrt(3)./3.*(2.*br+1))+ 1./2.*log(1+1./br+1./(br.*br))+sqrt(3).*pi./6;
    FF = 2.*rval.^2./T(xx,yy).*F;
    res = res + FF.*((1/2).*((-1).*utdt(xx,yy)+2.*ut(xx,yy).^2.*utdt(xx,yy)+(-1).*ut(xx,yy).^4.*utdt(xx,yy)+ut(xx,yy).*utdx(xx,yy).*ux(xx,yy)+(-1).*ut(xx,yy).^3.*utdx(xx,yy).*ux(xx,yy)+(-1).*ut(xx,yy).*ux(xx,yy).*uxdt(xx,yy)+ut(xx,yy).^3.*ux(xx,yy).*uxdt(xx,yy)+uxdx(xx,yy)+(-1).*ut(xx,yy).^2.*uxdx(xx,yy)+ux(xx,yy).^2.*uxdx(xx,yy)+ut(xx,yy).^2.*ux(xx,yy).^2.*uxdx(xx,yy)+ut(xx,yy).*utdy(xx,yy).*uy(xx,yy)+(-1).*ut(xx,yy).^3.*utdy(xx,yy).*uy(xx,yy)+ux(xx,yy).*uxdy(xx,yy).*uy(xx,yy)+ut(xx,yy).^2.*ux(xx,yy).*uxdy(xx,yy).*uy(xx,yy)+(-1).*ut(xx,yy).*uy(xx,yy).*uydt(xx,yy)+ut(xx,yy).^3.*uy(xx,yy).*uydt(xx,yy)+ux(xx,yy).*uy(xx,yy).*uydx(xx,yy)+ut(xx,yy).^2.*ux(xx,yy).*uy(xx,yy).*uydx(xx,yy)+uydy(xx,yy)+(-1).*ut(xx,yy).^2.*uydy(xx,yy)+uy(xx,yy).^2.*uydy(xx,yy)+ut(xx,yy).^2.*uy(xx,yy).^2.*uydy(xx,yy)));
    
  end

  function phorizon()
    figure(fighorizon);
    [xp,yp] = ndgrid(1:GRID_SIZE, 1:GRID_SIZE);
    size(xp)
    size(yp)
    size(horizon)
    surf(xp, yp, horizon);
    zlim( [5-4 5+4])
    caxis([5-4 5+4])
  end

  function pmetric2()
    figure(figmetric);
    plotsx = 2; 
    plotsy = 2;

    plotmag = 100;


    mdata = reshape(gr(GRID_SIZE/2, :,:), GRID_SIZE, R_SIZE);
    subplot(plotsy, plotsx, 1);
    image(mdata*CMAP_OFFSET/1000+CMAP_OFFSET);
    set(gca, 'XTick', rticks);
    set(gca, 'XTickLabel', rticklabels);

    subplot(plotsy, plotsx, 2);
    image(mdata*CMAP_OFFSET/100000+CMAP_OFFSET);
    set(gca, 'XTick', rticks);
    set(gca, 'XTickLabel', rticklabels);

    subplot(plotsy, plotsx, 3);
    image(mdata*CMAP_OFFSET/1000000+CMAP_OFFSET);
    set(gca, 'XTick', rticks);
    set(gca, 'XTickLabel', rticklabels);

    subplot(plotsy, plotsx, 4);
    image(mdata*CMAP_OFFSET/10000000+CMAP_OFFSET);
    set(gca, 'XTick', rticks);
    set(gca, 'XTickLabel', rticklabels);

    if (VIDEO) 
      frame = getframe(gcf);
      writeVideo(vt,frame);
    end

  end
 
  function pmetric()
    figure(figmetric);
    plotsx = 3; 
    plotsy = 4;

    plotmag = 100;

    mdata = reshape(g(GRID_SIZE/2, :,:,1,1), GRID_SIZE, R_SIZE);
    subplot(plotsy, plotsx, 1);
    image(mdata*CMAP_OFFSET/plotmag+CMAP_OFFSET);
    set(gca, 'XTick', rticks);
    set(gca, 'XTickLabel', rticklabels);
    title('gtt');

    mdata = reshape(g(GRID_SIZE/2, :,:,1,2), GRID_SIZE, R_SIZE);
    subplot(plotsy, plotsx, 2);
    image(mdata*CMAP_OFFSET/plotmag+CMAP_OFFSET);
    set(gca, 'XTick', rticks);
    set(gca, 'XTickLabel', rticklabels);
    title('gtx');

    mdata = reshape(g(GRID_SIZE/2, :,:,1,3), GRID_SIZE, R_SIZE);
    subplot(plotsy, plotsx, 3);
    image(mdata*CMAP_OFFSET/plotmag+CMAP_OFFSET);
    set(gca, 'XTick', rticks);
    set(gca, 'XTickLabel', rticklabels);
    title('gty');

    mdata = reshape(g(GRID_SIZE/2, :,:,2,2), GRID_SIZE, R_SIZE);
    subplot(plotsy, plotsx, 4);
    image(mdata*CMAP_OFFSET/plotmag+CMAP_OFFSET);
    set(gca, 'XTick', rticks);
    set(gca, 'XTickLabel', rticklabels);
    title('gxx');

    mdata = reshape(g(GRID_SIZE/2, :,:,2,3), GRID_SIZE, R_SIZE);
    subplot(plotsy, plotsx, 5);
    image(mdata*CMAP_OFFSET/plotmag+CMAP_OFFSET);
    set(gca, 'XTick', rticks);
    set(gca, 'XTickLabel', rticklabels);
    title('gxy');

    mdata = reshape(g(GRID_SIZE/2, :,:,3,3), GRID_SIZE, R_SIZE);
    subplot(plotsy, plotsx, 6);
    image(mdata*CMAP_OFFSET/plotmag+CMAP_OFFSET);
    set(gca, 'XTick', rticks);
    set(gca, 'XTickLabel', rticklabels);
    title('gyy');

    mdata = reshape(g(GRID_SIZE/2, :,:,1,4), GRID_SIZE, R_SIZE);
    subplot(plotsy, plotsx, 7);
    image(mdata*CMAP_OFFSET/plotmag+CMAP_OFFSET);
    set(gca, 'XTick', rticks);
    set(gca, 'XTickLabel', rticklabels);
    title('gtr');

    mdata = reshape(g(GRID_SIZE/2, :,:,2,4), GRID_SIZE, R_SIZE);
    subplot(plotsy, plotsx, 8);
    image(mdata*CMAP_OFFSET+CMAP_OFFSET);
    set(gca, 'XTick', rticks);
    set(gca, 'XTickLabel', rticklabels);
    title('gxr');

    mdata = reshape(g(GRID_SIZE/2, :,:,3,4), GRID_SIZE, R_SIZE);
    subplot(plotsy, plotsx, 9);
    image(mdata*CMAP_OFFSET+CMAP_OFFSET);
    set(gca, 'XTick', rticks);
    set(gca, 'XTickLabel', rticklabels);
    title('gyr');
    %if (GRID_SIZE >= 64)
    %  truesize(figmetric);
    %end
    mdata = reshape(detg(GRID_SIZE/2, :,:), GRID_SIZE, R_SIZE);
    subplot(plotsy, plotsx, 10);
    image(mdata*CMAP_OFFSET/1000+CMAP_OFFSET);
    set(gca, 'XTick', rticks);
    set(gca, 'XTickLabel', rticklabels);
    title('detg');

    mdata = reshape(gr(GRID_SIZE/2, :,:), GRID_SIZE, R_SIZE);
    subplot(plotsy, plotsx, 11);
    image(mdata*CMAP_OFFSET/100000+CMAP_OFFSET);
    set(gca, 'XTick', rticks);
    set(gca, 'XTickLabel', rticklabels);
    title('gr');

  end

  function r = phydro()
    frame_active = (mod(t-1,int16(1/Dt/SECONDS_PER_LIGHTYEAR/FPS)) == 0) || EVERY_FRAME;
    if VIDEO & ~frame_active 
      return
    end
    if VIDEO
      disp(t);
    end
    plotsx = 4;
    plotsy = 3;

    figure(fighydro);

    subplot(plotsy, plotsx, 1);
    image((T-TEMP_OFFSET)*CMAP_OFFSET +CMAP_OFFSET); 
    title('T');
    
    % if (VIDEO) & (frame_active)
    %   frame = getframe;
    %   writeVideo(vt,frame);
    % end

    subplot(plotsy, plotsx, 2);
    image(ux*CMAP_OFFSET +CMAP_OFFSET);
    title('ux');
    % if (VIDEO) & (frame_active)
    %   frame = getframe;
    %   writeVideo(vvx,frame);
    % end

    subplot(plotsy, plotsx, 3);
    image(uy*CMAP_OFFSET +CMAP_OFFSET);
    title('uy');
    % if (VIDEO) & (frame_active)
    %   frame = getframe;
    %   writeVideo(vvy,frame);
    % end

    subplot(plotsy, plotsx, 4);
    image(sqrt(ux.*ux + uy.*uy)*CMAP_OFFSET*10 +CMAP_OFFSET);
    title('sqrt(ux^2+uy^2)');
    % if (VIDEO) & (frame_active)
    %   frame = getframe;
    %   writeVideo(vv,frame);
    % end

    subplot(plotsy, plotsx, 5);
    image(Tdt+CMAP_OFFSET);
    title('dtT');

    subplot(plotsy, plotsx, 6);
    image(uxdt*CMAP_OFFSET/32+CMAP_OFFSET);
    title('dtux');
    subplot(plotsy, plotsx, 7);
    image(uydt*CMAP_OFFSET/32+CMAP_OFFSET);
    title('dtuy');

    subplot(plotsy, plotsx, 9);
    image(P1*CMAP_OFFSET/8+CMAP_OFFSET);
    title('dtP1');
    subplot(plotsy, plotsx, 10);
    image(P2*CMAP_OFFSET/8+CMAP_OFFSET);
    title('dtP2');


    title(['t = ',num2str(t*Dt)]);

    if (GRID_SIZE >= 64)
      truesize(fighydro);
    end

    %if (VIDEO) & (frame_active)
    %  frame = getframe(gcf);
    %  writeVideo(vt,frame);
    %end
  end

  % [-1 0 1
  function r = dx2(arg)
    r = zeros(GRID_SIZE, GRID_SIZE);
    for y = 1:GRID_SIZE
      r(1,y) = - arg(GRID_SIZE,y)/2 + arg(2, y)/2;
    end
    for y = 1:GRID_SIZE
      r(GRID_SIZE,y) =  -arg(GRID_SIZE-1,y)/2 + arg(1, y)/2;
    end
    for x = 2:GRID_SIZE-1
      for y = 1:GRID_SIZE
        r(x, y) = (arg(x+1,y) - arg(x-1,y))/2;   
      end
    end 
    r = r / DX;
  end
  function r = dy2(arg)
    r = zeros(GRID_SIZE, GRID_SIZE);
    for x = 1:GRID_SIZE
      r(x,1) = - arg(x,GRID_SIZE)/2 + arg(x,2)/2;
    end
    for x = 1:GRID_SIZE
      r(x,GRID_SIZE) =  -arg(x,GRID_SIZE-1)/2 + arg(x,1)/2;
    end
    for x = 1:GRID_SIZE
      for y = 2:GRID_SIZE-1
        r(x, y) = (arg(x,y+1) - arg(x,y-1))/2;   
      end
    end 
    r = r / DX;
  end
  function r = dxdx2(arg)
    r = zeros(GRID_SIZE, GRID_SIZE);
    for y = 1:GRID_SIZE
      r(1,y) = arg(GRID_SIZE,y) + arg(2, y) - 2*arg(1,y);
    end
    for y = 1:GRID_SIZE
      r(GRID_SIZE,y) =  arg(GRID_SIZE-1,y) + arg(1, y) - 2*arg(GRID_SIZE,y);
    end
    for x = 2:GRID_SIZE-1
      for y = 1:GRID_SIZE
        r(x, y) = arg(x+1,y) + arg(x-1,y) - 2*arg(x,y);
      end
    end 
    r = r / DX / DX;
  end
  function r = dydy2(arg)
    r = zeros(GRID_SIZE, GRID_SIZE);
    for x = 1:GRID_SIZE
      r(x,1) = arg(x,GRID_SIZE) + arg(x,2) - 2*arg(x,1);
    end
    for x = 1:GRID_SIZE
      r(x,GRID_SIZE) =  arg(x,GRID_SIZE-1) + arg(x,1) - 2*arg(x,GRID_SIZE);
    end
    for x = 1:GRID_SIZE
      for y = 2:GRID_SIZE-1
        r(x, y) = arg(x,y+1) + arg(x,y-1) - 2*arg(x,y);   
      end
    end 
    r = r / DX / DX;
  end

end


  function r = dx(arg)
    global d1;
    r = d1 * arg;
  end
  function r = dy(arg)
    global d1;
    r = arg * (d1.');
  end
  function r = dxdx(arg)
    global d2;
    r = d2 * arg;
  end
  function r = dydy(arg)
    global d2;
    r = arg * (d2.');
  end
